function [sj_out, p_out, tar_s] = new_eq_sj(p_ini,mc_cf,theta2,ns,ind_market,x2,tax_ava,s_tar, delta, A1, QI,year,FC,total_market,ss)
sel_mrkt=ind_market(year>=2019);
% 
options = optimset('Algorithm','sqp','GradObj', ...
'off','GradCon','off','DerivativeCheck','off', 'maxiter', 4000,...
'MaxFunEvals', 1000000,'TolFun',1e-12,'TolX',1e-12);


%p_out=p_ini/1000;
p_out=p_ini;
tar_s=ones(max(sel_mrkt),1).*0.35;
p_2={};
tx={};
v_c=getranddraw(length(theta2)-1, ns);


expmu1=expmu(theta2,v_c,x2,p_out);
sj_out=mktshares(delta,expmu1,QI,ind_market);

parfor ii=min(sel_mrkt):max(sel_mrkt)
    ii
    temp=find(ind_market==ii);
    ind_market_i=ind_market(temp);
    [~,tax_ava_i]=ismember(tax_ava,temp);
    tax_ava_i(tax_ava_i==0)=[];
    p_nt_i=p_ini(temp);
    FC_i=FC(temp,ss);
    delta_i=delta(temp);
    mc_i=mc_cf(temp);
    A1_i=A1(temp,temp);
    x2_i=x2(temp,:);
    ini_var=[s_tar;p_nt_i];
    mrks_i=repmat(total_market(ii),length(temp),1);
    QI_i=QI(temp,:);
    
    lb=[0; 0*ones(length(ini_var)-1,1)];%lower bound for price and import tariff
    ub=[0.19; inf*ones(length(ini_var)-1,1)]; %upper bound for proce and import tariff
    
    tic
    [out_1,~]=fmincon(@(ini)new_sj(ini,delta_i,x2_i,theta2,v_c,tax_ava_i,FC_i,mc_i,A1_i,mrks_i),ini_var,[], [], [], [], ...
    [lb],[ub],[], options);
    toc
    
    cf_salestax=ones(size(delta_i))*0.19;
    tariff_cf=(ones(size(delta_i))*0.35);
    tariff_cf(tax_ava_i)=out_1(1);
    mc_cc=mc_i.*(1+tariff_cf);
    expmu1=expmu(theta2,v_c,x2_i,out_1(2:end).*(1+cf_salestax));
    s_indiv=ind_shares(exp(delta_i),expmu1); 
    sj_cf=mean(s_indiv,2);
    price1=(mc_cc+markup_sim(theta2,v_c,delta_i,s_indiv,A1_i,sj_cf)).*(1+cf_salestax);
    expmu2 = expmu(theta2,v_c,x2_i,price1);
    sj_2=mktshares(delta_i,expmu2,QI_i,ind_market_i);
    
    s_2{ii}=sj_2
    p_2{ii}=price1;
    tx{ii}=out_1(1);
    
end

for jj=min(sel_mrkt):max(sel_mrkt)
    temp2=find(ind_market==jj);
    p_out(temp2)=p_2{jj};
    tar_s(jj)=tx{jj};
    sj_out(temp2)=s_2{jj};
end


end% function computing counterfactual price equilibrium  

